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In this paper we present the implementation of an efficient formalism for the generation of arbitrary 
non-Gaussian initial conditions for use in N-body simulations. The methodology involves the use of 
a separable modal approach for decomposing a primordial bispectrum or trispectrum. This approach 
allows for the far more efficient generation of the non-Gaussian initial conditions already described 
in the literature, as well as the generation for the first time of non-separable bispectra and the 
special class of diagonal-free trispectra. The modal approach also allows for the reconstruction of 
the spectra from given realisations, a fact which is exploited to provide an accurate consistency 
check of the simulations. 



INTRODUCTION 

Testing for deviations of primordial density fluctua- 
tions from Gaussianity represents one of the most ac- 
tive areas of research in cosmology today (see for exam- 
ple [HUT]). Detection of an appreciable deviation would 
violate the current slow-roll inflationary paradigm. To 
date most tests of non-Gaussianity have focussed on con- 
straining the primordial skewness, described by the three- 
point function or bispectrum, using the cosmic microwave 
background (CMB). The resultant CMB non-Gaussianity 
may be simply related to its primordial 'seed' via trans- 
fer functions. This relationship reflects the fact that the 
CMB is well described by linear theory. Large scale struc- 
ture (LSS) as a three dimensional data source, unlike the 
two dimensional CMB, offers the possibility of a vast im- 
provement in constraining non-Gaussianity. However, a 
major drawback is the non-linear relation between the 
primordial density fluctuation and the resulting distri- 
bution of structure. For this reason, the investigation 
of non-Gaussianity using LSS must take a more empir- 
ical approach relying on N-body simulations. Owing to 
the complexity involved in generating non-Gaussian ini- 
tial conditions, relatively few models have been studied 
to date. In fact, aside from the local model only the 
non-Gaussian bispectra of the equilateral and orthogonal 
shapes have been studied [HI [J3] . The implementation 
in these latter cases involved an extremely computation- 
ally expensive algorithm. In this paper we describe an 
efficient method to create non-Gaussian initial conditions 
for arbitrary bispectra and the special class of diagonal- 
free trispectra. The approach makes use of the separable 
decomposition of the primordial spectra, which has been 
exploited to considerable success in the case of the CMB 
[14H20j . In this paper we present a brief overview of the 
formalism (for a more detailed exposition see ref. [21]). 
We detail a non-trivial check of the simulations, verifying 
the accuracy and consistency of the approach. Finally we 
summarise our findings. 



ALGORITHM 

In this section we describe briefly the algorithm for the 
generation of non-Gaussian initial conditions. We assume 
that the density field is statistically isotropic. Our treat- 
ment is universal in that it covers general bispectra and 
the class of trispectrum models which depend only on the 
magnitude of its wavenumbers, i.e. diagonal-free trispec- 
tra. This case covers almost all trispectra discussed to 
date in the literature, except for the diagonal-dependent 
local (tnl) trispectrum. However, the local tjvx may 
be simply generated using the following expansion about 
two Gaussian fields, 4>g and tpa, (where (^g^g) = 
[22]) 1 

(=<f> G + lP G + f NL (<(>%- (4» G )). (1) 

The algorithm described here incorporates the generation 
of an explicit trispectrum in the absence of a bispectrum 
and vice versa. It should be noted that the bispectrum 
term also generates an implicit trispectrum. In the case 
of the local model this 'spurious' trispectrum is the t/vl 
model described above. Such trispectra may not be de- 
sirable in other models and so should be subtracted out 
[21] , This issue will be addressed further in a future pa- 
per. 

Bispectrum 

As described in (21] . an arbitrary primordial bispec- 
trum, B{ki,k 2 ,k^), may be simulated by evaluating the 
quantity 

$B(k)= / ^^^(k-k'-k")$ G (k')<i> G (k") 

x B{k,k',k") 

P{k)P(k') + P{k)P[k") + P{k')P{k"Y 1 ' 



In the single field case ipQ is set to zero and tjvl = (6/jvz,/5) 2 . 
In general tjvi, obeys the inequality tjvi, > (6/jvi/5) 2 /2 |23| . 
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where $ G is a Gaussian random field with the required 
power spectrum P{k). This expression, written in con- 
volved form, was used to tackle some specific separable 
bispectrum models in refs. [T21 [T3] . It is directly related 
to that employed for creating non-Gaussian CMB map 
simulations [4 which was generalised with modal meth- 
ods in ref. |16j . The modal approach eliminated poten- 
tial non-Gaussian contributions to the CMB power spec- 
trum. Here, however, the power spectra in the denomina- 
tor must also be symmetrised to mitigate against these 
spurious effects [15] , The expression ^ is the natural 
choice for initial conditions since, for the local model of 
inflation, this procedure reduces to the usual convolution 
<j) G * <J) G . The primordial perturbation, $, given by 



(3) 



then obeys (in the limit of weak non-Gaussianity) the 
desired relations 

($(k 1 )$(k 2 )) = (27T) 3 S D (k 1 +k 2 )P(fcl), 

($(k 1 )$(k 2 )$(k 3 )) = (27r) 3 5 D (Y / K)F NL B(k 1 ,k 2 ,k 3 ). 

(4) 

The direct calculation of initial conditions via this pre- 
scription is not efficient in general due to the non- 
separable form of the integrand on the second line of 
([2]). However, this term may be rewritten in a separable 
form using the modal techniques described in |14H16j . In 
particular, we may expand the integrand within ([2| in 
the form 



B(k, k', k" 



P(k)P(k') + P(k)P(k") + P(k')P(k") 
J2^ s Mk)q s (k')q t (k" 



(5) 



needed to accurately account for most of the models de- 
scribed in the literature. We also express the Dirac delta 
function in the form 



Mk) 



1 



(2tt) e 



d 3 xe lkx . 



(6) 



The bispectrum contribution may now be efficiently eval- 
uated as 



$ B (k) 



x%q {r (k) 



<M s (x)M t} (x), (7) 



where the filtered density perturbations, M s (k), are 
given by 



M,(x) 



d 3 k 

(2tt) 3 



$ G (k)g s (fc)e- 



(8) 



Thus the evaluation has been reduced to the calculation 
of a series of fast Fourier transforms. These expressions 
are to be evaluated in a box corresponding to a maximum 
wavenumber fc max - Care must be taken to account for 
unwanted realisations of the discretisation of the Dirac 
delta function when the wavevectors, k^, align. This can 
be accounted for simply by restricting the range of the 
wavevectors to |ki| < 2/c max /3 for the calculation of $ B . 
This limitation is more than offset by the vast improve- 
ment in numerical speed and accuracy that the modal 
method offers. 

Once the non-Gaussian primordial potential <f> s (k) is 
generated it can be translated into the linear density per- 
turbation <5k, 2 at some initial redshift z using the Poisson 
equation and transfer function T{k). From 5^. z one can 
get initial particle positions and velocities for N-body 
codes using the Zel'dovich approximation [24] or second- 
order Lagrangian perturbation theory [35] 12S] ■ 



where the q r are one dimensional orthogonal polynomi- 
als on the domain of validity of the bispectrum, that is, 
the tetrahedral region prescribed by the closure condition 
imposed by the Dirac delta function. Note that the form 
of these mode functions q r is not important - whether 
polynomial, trigonometric, wavelet, etc - provided they 
form a complete set (those used in this paper are close to 
Lcgcndrc polynomials and arc defined in ref. 16J). We 
may introduce a partial ordering on the indices used in 
their 3D products and write ^2 rst af st q r (k)q s (k')qt(k") = 
Yjn={prs} a nq{r{k)q s {k')qt}{k"), where {...} represents 

the symmetrised quantity 2 . The coefficients a® charac- 
terise the specific model under scrutiny. As has been 
shown in [IS], relatively few modes (n max = O(30)) are 



2 In what follows we use the compact notation Q n (k,k' ,fc") to 
represent q^ r (k)q a (k')q t y(k"). 



Trispectrum 

As indicated already, we shall only consider the spe- 
cial class of diagonal-free trispectra in this paper. Such 
trispectra are given by the following four-point connected 
correlator 

($(k 1 )$(k 2 )$(ki)$(ki)) c =(2^) 3 Mk! + k 2 + k 3 + k 4 ) 

x G NL T(k 1 ,k 2 ,k 3 ,k 4 ). 

(9) 

A primordial perturbation with the correct power spec- 
trum and trispectrum is then given by 



(10) 



where 



$T(k) = 1 dk ^ )6 d k S D (k k> k" - k'") 

T(k, k', k", k'") ^G^G^^G^uy 



P(k)P{k')P(k") + 3 perms 



(11) 
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Without the use of separable methods this integral would 
be intractable in general. However, separable methods 
outlined in [TTl [19] may again be used to greatly simplify 
the calculation. In particular, we may write 



T i^k ^ k . k . k ^ 



P(k)P(k')P{k") + 3 perms 



(12) 



where the q r are one dimensional orthogonal polynomials 
on the domain of validity of the trispectrum, and where 
m represents the partial ordering m = {rstu}. In what 
follows we may refer to the quantity Q m to represent 
Q{r<jsQt<iu} -The trispectrum contribution now becomes 



<f T (k) = 



q {r (k) 



<M s (x)M t (x)M u} (x), 
(13) 



where the filtered perturbations, M s (x), are as in equa- 
tion Q except for the replacement of q s by q s . Avoid- 
ance of unwanted images of the Dirac delta restricts the 
domain of validity to k^ < /c max /2. 



ALGORITHM VALIDATION 

In order to test the accuracy of the algorithm employed 
it is necessary to establish the convergence of the aver- 
age of estimators of particular realisations to the expec- 
tation value of the estimator. It should be noted that 
the efficacy of the primordial decomposition has been 
tested thoroughly in [TH [TSl [T^ with accuracy of at least 
O(90 — 95%) achievable in the case of the bispectrum 
with < 0(30) modes and in the case of the trispectrum 
with < 0(50) modes. 



Bispectrum estimation 

An estimator for the bispectrum is given by |21j 
II? =1 d 3 k 4 fo(kj + k 2 + k 3 )B(h,k 2) k 3 ) 



£ = 



(2tt) 6 P(Ai)P(fe)P(fe 3 ) 

x [0 kl 0k 2 0k 3 -3(0 kl k2 )0 k3 



(14) 



The expectation value of this estimator is given by 



v 



dk\dk 2 dk; 



V B 



kik 2 k 3 B 2 (ki,k 2 , k 3 
P(h)P(k 2 )P(k 3 ) 



-, (15) 



where Vb is the tetrahedral domain allowed by the trian- 
gle condition on the wavenumbers ki, and V is a volume 



factor given by V = (27t) 3 <5d(0). Again we expand the 
theoretical bispectrum in a separable form 



y/kik 2 k 3 B(ki,k 2 ,k 3 ) Q' ru ^ ru \ ru \ 

(16) 

This quantity is different to ^ used for the initial con- 
ditions because here it represents an expansion of the 
predicted signal-to-noise for the given bispectrum model. 
However, the two sets of expansion coefficients a® and 
a® can be directly related for any bispectrum using a 
matrix transformation, so we only need calculate one set 



of coefficients. The expansion ( 16 1 allows us to write the 



estimator and its expectation value in the form |21j 

£ =Y. a n /rf 3 x[Mr(x)A/ s (x)M t (x) 

n 

-<M {r (x)M.(x))M 1} (x)], (17) 
(£)=5>fa£' 7nm , (18) 



where 

r Ynm 

M r (x) 



V 



* -'Vb 
3 



dk\dk 2 dk 3 Q n (k\,k 2 , k 3 )Q m (ki,k 2 , fc 3 ), 



cPk kgr (fc) ik> 
(2tt)3 ^/kPjk) 6 



Expressing £ = J2 n a n fin > wu -h fin defined by equa- 
tion (17), we establish that 



(fin 



E 



(19) 



It is convenient to create an orthonormal set of 
mode functions lZ n from the product functions Q n 
(this may be done using a Gram-Schmidt orthog- 
onalisation using the the inner product (fg) — 
J VB dk 1 dk 2 dk 3 f(k 1 ,k 2 ,k 3 )g(k 1 ,k 2l k 3 )). In terms of 
these mode functions the consistency relationship may 
be easily shown to give 



(20) 
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Trispectrum estimation 



Hence, we establish that 



In the case of the diagonal-free trispectrum, the esti- 
mator and its expectation value take the form |21j 



Wn) = E ^,mlnm <=► (Pn) =*f,m, < ™ > 



£ = 



ntirf 3 k t jpQci + k 2 + k 3 + k 4 )T(fci, fc 2 , fc 3 , fc 4 ) 
(2tt) 9 P(fc 1 )P(fc 2 )P(fc 3 )P(fc 4 ) 
x [0 kl 0k 2 0k 3 0k 4 -6(0 kl k2 )0 k3 k4 

+ 3(0 kl k2 )(0 k3 k4 )], (21) 

>M-|fc 2 3|), (22) 



^E fc »- l^34| 



where fc 34 = ki + k 2 ~ k 3 — k 4 and Vt represent the 
domain of validity of the wavenumbers fc; as imposed by 
the Dirac delta function. Testing for the accuracy of the 
initial conditions in this case is slightly more involved 
than the bispectrum test. In particular, we can achieve 
this by using separable expansions of the following two 
theoretical quantities, 



y/kik 2 k 3 k 4: T(ki, k 2 ,k 3 ,k 4 ) 
v /P(fc 1 )P(fc 2 )P(fc 3 )P(fc 4 ) 



Vfclfc2fc 3 fc4T(fei, fc 2 , fc 3 , fe 4 ) 

^P{ki)P{k 2 )P(k 3 )P{k 4 ) 



= E* 



M4 



) n (fci,fc 2 ,fc 3 ,fc 4 ). 

(24) 



The estimator and its expectation value may now be ex- 
pressed in the form 



£ = E^ Q 

(£) = ^ 1 a ?,n a 2,rn'Jnrni 



(25) 
(26) 



where 



13 Q 



d d x |M r (x) M s (x) M t (x) M u (x) 
-6(M {r (x)M s (x))M t (x)M u} 
+ 3 ( M {r (x) M s (x) ) (M t (x) M u} ) 



(27) 



V 



(2tt) 6 



Ilf =1 dkiQ n (ki,k 2 , k 3 , fc 4 ) 



M r (x) 



x Q m (fci,fc 2 ,/c 3 ,fc 4 ), (28) 
d 3 k k g T .(fc) pik . x 
(2tt) 3 7fcP(fc) e 



(29) 



where the superscript 7?. refers to coefficients with respect 
to the orthonormal mode functions lZ n created from the 
product functions Q n . 



RESULTS 

In order to demonstrate the efficacy of these modal 
methods we have generated non-Gaussian initial con- 
ditions for the following bispectrum models (see e.g. 
ref. [H]): local, equilateral, constant, orthogonal and 
flattened (non-separable case). In addition, we have 
created trispectrum initial conditions for the local g^L 
model and the equilateral (c x ) model, as well as the con- 
stant model (see ref. [EI]). The decomposition of the pri- 
mordial shapes, as described by equations ^ and ( 12 ) 
respectively, was calculated first. Using these expansion 
coefficients the initial conditions were generated. The ac- 
curacy of these initial conditions was then tested using 
the bispectrum and trispectrum estimation techniques 
described in the previous section. The primordial de- 
compositions against which the initial conditions were 
compared are described in the case of the bispectrum by 



eqn (16), and in the case of the trispectrum by eqn (24 1 




k [h/MPc] 

Figure 1: Plot of input power P(k) for $ G (red) and mea- 
sured powers of /jvl ( 1 >s /2 for the local shape with /jvl = 100 
(black), /ati,'E )S /2 for the equilateral shape with Jnl~ 400 
(blue) and $ G + f NL ® B /2 for the two cases (green). $ was 
calculated for 100 realisations of "3? G using Q on a 256 3 grid 
in a (100Mpc/h) 3 box. 
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Bispectrum Results 



The primordial decompositions ^ and ( 16 ) are eval- 
uated using n max = 30 modes in the case of the lo- 
cal, equilateral, constant and orthogonal models, while 
n mnx = 80 modes were required in the case of the (non- 
separable) flattened model in order to achieve a more 
accurate fit. Since the initial conditions are calculated 
in the absence of inhomogeneities we may evaluate the 
coefficients of the estimator, j3® , using only the term 
/ rf 3 xM r (x)M„(x)M,(x) in equation ( flf) ,. Evaluation of 
each initial condition simulation is an extremely efficient 
operation, with each M r calculated in parallel. In the 
case of a f024 3 grid, a full initial condition simulation 
can be generated in approximately 0(1) hour of compu- 
tational time using only 6 cores. However, the main pur- 
pose of this paper is to provide a proof of concept, demon- 
strating a practical implementation of the methodology 
presented, so simulations have been carried out using a 
smaller 256 3 grid, unless otherwise stated. 

Before presenting the bispectrum validation results, we 
note that the non-Gaussian contribution /^ ri ($ B $ s )/4 
to the power spectrum is small compared to the Gaussian 
power, as can be seen in Figure[T]for both local and equi- 
lateral models. This is important because the simulated 
non-Gaussian contribution in any prescription must not 
modify the underlying power spectrum. 

In Figure [2] we plot a comparison of the theoretical bis- 
pectrum modes ct% and the estimated bispectrum modes 
(0^) from fOO simulation realizations (la error bars are 
shown). The agreement between the theoretical predic- 
tion and the averaged simulations is striking. In order to 
establish the accuracy of the approach we present in Ta- 
ble [I] the correlations between the primordial shape and 
the primordial decompositions ([5| and (16 1, as well as 
the correlation between the decomposition ( 16 ) and the 



average of the realistions f3® . The amplitude of the bis- 
pectrum -F/vl is a l so given in the table in each case. It is 
clear from the table that given a particular decomposition 
the average of the realisations is almost exact. The only 
limitation is the number of modes chosen to perform the 
primordial decompositions. However, in each of the cases 
considered here, the accuracy of the theoretical decom- 
position is greater than 90%. It is important to note that 
the only limitations here have emerged in the theoretical 
domain, rather than in the simulated initial conditions 
which faithfully represent the decomposed bispectra they 
are given by (to better than 99%). This theory limita- 
tion is easily circumvented by extending out to further 
coefficients in the modal expansion or by adapting the 
underlying modes for the case under investigation. 



Model 


Yl 

' ''max 




Shape vs 
Decomp 
Eqn (B) 


Shape vs 
Decomp 
Eqn ( 16 \ 


(B Q ) vs 
Decomp 
Eqn (161 


Local 


30 


100 


100% 


92.6% 


99.5% 


Equil 


30 


200 


99.7% 


99.7% 


99.8% 


Const 


30 


200 


99.9% 


100% 


99.3% 


Orthog 


30 


200 


98.7% 


98.9% 


99.7% 


Flat 


80 


200 


91.8% 


90.6% 


99.1% 



Table I: Correlation between the primordial bispectrum shape 
and the modal decompositions |5| and ( |16[ ), as well as the 
correlation between the average of the realisations {fin) and 
the decomposition (16 1. 
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Figure 2: Plo t of the theoretical local bispectrum modes as 
described by ( 16 1 compared to the average of the estimated 
modes f}% from 100 realisations. The modes are compared in 
the rotated (orthonormal) 1Z basis. 



Trispectrum Results 

The primordial models considered in the case of the 
trispectrum, i.e. the local (gNh), equilateral (c\) and 
constant models, are quite well-behaved and can be ac- 



curately decomposed, as described by equations (12 1 and 
(24), using just n max = 18 modes. In order to calculate 

the coefficients f3 n it is necessary to calculate correla- 
tors of the form (M r (x)M s (x)) , as described by equation 
(27 1. Here, we have carried out fOOO simulations in order 



to measure this correlation accurately. This operation - 
to test the accuracy of the initial conditions - represents 
the most numerically intensive operation taking approx- 
imately 8 hours on 28 cores. Again a grid size of 256 3 
is used. The are then calculated and the average 
over 200 simulations is calculated and compared to the 
theoretical prediction (24). In Figure [3] the modes are 

compared, in the case of the local (gNh) model, with la 

—iz 

error bars included for the average of the f3 n . Clearly, 
the two sets of modes are highly correlated, validating 
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the trispectrum methodology. In Table |TT] we present 
the correlations between the primordial shape and the 



primordial decompositions (12) and (24), as well as the 



correlation between the decomposition (24) and the av- 
erage of the realisations /3®, with the amplitude of the 
trispectrum given by Gnl (see equation (10)). The al- 



most 100% correlation in each case verifies the accuracy 
and validity of the approach. It should be noted that 
while the choice of Gnl is arbitrary, a lower amplitude 
will, of course, result in a lower signal to noise for the 
estimator. 



Model 


^max 


Gnl 


Shape vs 
Decomp 
Eqn |[12| 


Shape vs 
Decomp 
Eqn p4| 


(P Q n ) vs 
Decomp 
Eqn (24| 


Local 


18 


5 x 10 6 


99.6% 


100% 


100% 


Equil 


18 


5 x 10 6 


99.4% 


99.8% 


100% 


Const 


18 


5 x 10 6 


99.9% 


99.9% 


100% 



Table II: Correlation bet wee n the primordial shape and the 
decompositions (|5| and (16 1, as well as the correlation be- 
tween the average of the realisations (Pn) and the decompo- 
sition JT6l). 



models, as well as the (non-separable) flattened model. 
We have also presented the implementation for the class 
of diagonal- free trispectra, including equilateral (ci), lo- 
cal (<7atl) and constant models. If primordial bispectra 
and trispectra are to be simulated together, then it may 
be necessary for a spurious primordial bispectrum contri- 
bution to the trispectrum to be subtracted out (except in 
the simplest local tnl model). However, the subtraction 
of such terms involves extending the work presented here 
to 'diagonal-dependent' trispectra. We defer a detailed 
quantitative analysis of this trispectrum issue to a future 
publication. Nonetheless, the present work represents a 
significant step forward in opening up the efficient inves- 
tigation of primordial non-Gaussianity using large scale 
structure. 

For the non-Gaussian initial conditions described in 
this paper, a non-trivial consistency check has been car- 
ried out to verify their efficacy. It has been established 
that - once unwanted images of the Dirac delta function 
are accounted for - the algorithm employed here is unbi- 
ased and accurate. 
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Figure 3: Plot of the theoretical modes as described by (24) 

compared to the average of the modes, over 200 realisa- 
tions for the local g^L model. The modes are compared in 
the rotated (orthonormal) TZ basis. 



SUMMARY 

We have described the generation of non-Gaussian ini- 
tial conditions for use in N-body simulations. Exploiting 
the use of a separable modal approach to analyse primor- 
dial non-Gaussianity, the algorithm is reduced to a series 
of fast Fourier transforms and a three dimensional inte- 
gral. We have described an application of the approach to 
arbitrary bispectra, presenting for brevity the implemen- 
tation of the local, equilateral, constant and orthogonal 
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